
restore,'./local/BGPS_V1_pvec.sav',/ver
restore,'./local/BGPS_V1_emaf.sav',/ver

eind = where(constrain.emaf)
pvec = pvec[eind]
constrain = constrain[eind]

near = pvec.stat.kdar EQ 'N'
far  = pvec.stat.kdar EQ 'F'
tan  = pvec.stat.kdar EQ 'T'
print,'NEAR: ',fix(total(near)),' FAR: ',fix(total(far)),$
      ' TAN: ',fix(total(tan))

;; Only include well-constrained!!!1!
wcind = where(constrain.post, nwc)
print,'NWC: ',nwc
duse  = pvec[wcind].stat.duse / 1.d3
emaf = emaf[wcind]

mw = omni_read_conffile('./conffiles/galactic_params.conf')

print,'DUSE: ',size(duse,/dim)

full_duse  = duse
full_emaf = emaf

myps,'./emaf_paper/plots/distance_comparison_paper_f13a_bw.eps',$
     xsize=3.1,ysize=3.1,/bw

;; multiplot,[2,1],xgap=0.035

symper = 0.8
cs = 0.7

masercolor = 'blk4';'Magenta'
nirexcolor = 'blk6';'Cyan'

range=[0,15]

Plot,[0],[0],/isotropic,charsize=cs,/nodata,xr=range,yr=range,/xst,/yst,$
     xtit='Derived Heliocentric Distance  [kpc]',$
     ytit='Comparison Set Distance  [kpc]',xmargin=[8,1],ymargin=[4,1]


;; Plot the 'alley'
cgPlots,0,1
cgPlots,range[1]-1-0,range[1]-0,/cont,linestyle=2,color='BLK5'
cgPlots,1,0
cgPlots,range[1]-0,range[1]-1-0,/cont,linestyle=2,color='BLK5'


;;=============================================================================
;;=============================================================================
;;=============================================================================
;; Do the GRS stuff -- create GRS structure
;; Load in necessary files
restore,'./ancillary/grs_kda_data.sav',/ver
n = nwc
angsep = 10.
emafgrs = replicate( {has:0b,$
                      clemdist:0.,$
                      v:0.,$
                      l:0,$
                      b:0.,$
                      cm21:'',$
                      dist:0.}, n)

;; Create arrays used for GRS identification
dsep = fltarr(n)
vsep = fltarr(n)
grsi = lonarr(n)

;; Loop through EMAF to check for GRS co-incidence
message,'Checking for GRS co-incidence...',/inf
FOR i=0L, n-1 DO BEGIN
   gcirc, 2, emaf[i].l, emaf[i].b, grs.glon, grs.glat, dis
   mind = min(dis,gind)
   vsep[i] = abs(emaf[i].vlsr - grs[gind].vlsr)
   
   ;; A DSEP =-1 means there is no coincidence for this object
   dsep[i] = (vsep[i] LE grs[gind].delv*1.0) ? mind/60. : -1
   IF dsep[i] GT angsep THEN dsep[i] = -1
   grsi[i] = gind
   
ENDFOR

emafgrs.has = dsep GE 0.
message,'There are '+string(total(emafgrs.has),format="(I0)")+$
                         ' GRS-associated sources.',/inf

;; Load emaf
emafgrs.clemdist = grs[grsi].dis
emafgrs.v        = grs[grsi].vlsr
emafgrs.l        = grs[grsi].glon
emafgrs.b        = grs[grsi].glat
emafgrs.cm21     = grs[grsi]._21cm_

;; Derive GRS distance based on Reid (2009) rotation curve

;; Read in GRS distance correction file
readcol,'./emaf_paper/corrected_grs_distance.txt',bad_gd_cnum,new_gd,format='I,F',$
        COUNT=ncorrected
IF ncorrected GT 0 THEN BEGIN
   ;; Sort GRS list
   gind        = SORT(bad_gd_cnum)
   bad_gd_cnum = bad_gd_cnum[gind]
   new_gd      = new_gd[gind]
   ;; Figure out where, etc...
   e_ind = WHERE_ARRAY(bad_gd_cnum,emaf.cnum,n_e)
   c_ind = WHERE_ARRAY(emaf.cnum,bad_gd_cnum,n_c)
   IF n_e NE n_c THEN message,'Drink some of that stout, chubby!'
   ;; IF n_e NE ncorrected THEN message,'Drink some of that porter, chubby!'
   IF n_e NE 0 THEN emafgrs[e_ind].clemdist = new_gd[c_ind]
ENDIF

;; Calculate tangent distances & how many at FAR KDIST
grs_dtan = MW.R0 * cos(emafgrs.l * !dtor) / cos(emafgrs.b * !dtor)/1.d3

;; If within 0.25 kpc of dtan, assign to dtan
TAN_BUFF = 0.0
grs_tan = abs(emafgrs.clemdist-grs_dtan) LT TAN_BUFF AND emafgrs.has
grs_far = (emafgrs.clemdist GE grs_dtan+TAN_BUFF) AND emafgrs.has
message,'N GRS at TANGENT d: '+string(long(total(grs_tan)),format="(I0)"),/inf
message,'N GRS at far KDIST: '+string(long(total(grs_far)),format="(I0)"),/inf

;; Calculate R09 distance based on vlsr & grs_far
FOR i=0L,n-1 DO BEGIN
   
   ;; Calculate grs_dist based on GRS KDA resolution
   IF emafgrs[i].has THEN $
      IF ~ grs_tan[i] THEN $
         emafgrs[i].dist = (grs_far[i]) ? $
                            KDIST(emafgrs[i].l,emafgrs[i].b,emafgrs[i].v,$
                                  /FAR)/1.d3 : $
                            KDIST(emafgrs[i].l,emafgrs[i].b,emafgrs[i].v,$
                                  /NEAR)/1.d3 $
      ELSE emafgrs[i].dist = MW.R0 * cos(emaf[i].l * !dtor) / $
                              cos(emaf[i].b * !dtor)/1.d3
ENDFOR

;;=============================================================================
;;=============================================================================
;;=============================================================================

print,'STATS: '
print,m4_stat(emafgrs.dist)
print,m4_stat(emafgrs.clemdist)
print,m4_stat(duse[0,*])

gind = where(emafgrs.has, ngrs)
print,'NGRS: ',ngrs

duse = duse[*,gind]
emaf = emaf[gind]
emafgrs = emafgrs[gind]

;; Find mismatches:
bind = where( abs(duse[0,*]-emafgrs.dist) GT 1.0, nb)
print,'NBAD: ',nb,' Matching rate = '+$
      string(float(ngrs-nb)/ngrs*100.,format="(F0.1,'%')")

;; Look at contrast limits
FOR climit = 0.1,0.31,0.05 DO BEGIN
   icont = where(emaf.c GE climit, ncont)
   bci = where( abs(duse[0,*]-emafgrs.dist) GT 1.0 AND emaf.c GE climit, nbc)
   print,'Contrast C > '+string(climit,format="(F0.2)"),$
         '  NBAD: ',nbc,' Matching rate = '+$
         string(float(ncont-nbc)/ncont*100.,format="(F0.1,'%')")
ENDFOR

;; Write out CNUMS for examination
openw,lun,'./emaf_paper/exam_grs_mismatch.txt',/GET_LUN
for ll=0,nb-1 do printf,lun,emaf[bind[ll]].cnum
close,lun
free_lun,lun

oploterror,duse[0,bind],emafgrs[bind].dist,duse[1,bind],duse[1,bind]*0.,/lobar,$
           psym=3,errcolor='BLK5'
oploterror,duse[0,bind],emafgrs[bind].dist,duse[2,bind],duse[2,bind]*0.,/hibar,$
           psym=3,errcolor='BLK5'

cgOplot,duse[0,*],emafgrs.dist,psym=16,symsize=0.35

ilo = where((duse[0,*]-emafgrs.dist) GT 1.0,nlo)
print,'BR: '
FOR iii=0,nlo-1 DO $
   print,emaf[ilo[iii]].cnum,duse[0,ilo[iii]],emafgrs[ilo[iii]].dist,$
         emaf[ilo[iii]].vlsr,emafgrs[ilo[iii]].v
ihi = where((duse[0,*]-emafgrs.dist) LT -1.0,nhi)
print,'UL: '
FOR iii=0,nhi-1 DO $
   print,emaf[ihi[iii]].cnum,duse[0,ihi[iii]],emafgrs[ihi[iii]].dist,$
         emaf[ihi[iii]].vlsr,emafgrs[ihi[iii]].v


;; inh3 = where(emaf.v_flag EQ 8)
;; cgOplot,duse[inh3,0],emaf[inh3].grs_dist,psym=15,color='brown'
;; ib4 = where( abs(emaf.b) GE 0.4)
;; cgOplot,duse[ib4,0],emaf[ib4].grs_dist,psym=15,color='brown'

;; ;; Check on 21-cm continuum & changed distance
;; i21 = where(emafgrs.21cm EQ strlowcase(strtrim('y',2)), n21)
;; print,'N21: ',n21
;; cgOplot,duse[i21,0],emaf[i21].grs_dist,psym=15,symsize=1.0,color='brown'

;; readcol,'./emaf_paper/corrected_grs_distance.txt',cnum,d,format='i,f'
;; e21 = emaf[i21]
;; ich = where_array(cnum,e21.cnum,nch)
;; print,'N changed when they should not have been: ',nch

testc = [3175,3355,3542,3738,4223,4338,4539,4544,5005,5237,5238,5289,5767,5868,6200,6277]

cin = where_array(testc,emaf.cnum,ncin)
print,'NCIN: ',ncin

print,'============================'
print,'Contrast Stats: ',m4_stat(emaf[cin].c)


duse = full_duse
emaf = full_emaf


;;=======================================
;; PLOT MARSHALL ET AL. (2009) NIREX DISTANCES
m09 = read_mrt('./ancillary/Marshall2009_Table1.txt')

m09ind = lonarr(nwc)
m09dist = fltarr(nwc)
FOR i=0L, nwc-1 DO BEGIN
   gcirc,2,emaf[i].l,emaf[i].b,m09.glon,m09.glat,dis
   d_min = min(dis,mind)
   m09ind[i] = mind
   
   IF d_min LE 1.*60. THEN BEGIN
      print,d_min, duse[0,i],m09[mind].dist,m09[mind].dist-duse[0,i]
      m09dist[i] = m09[mind].dist
   ENDIF
ENDFOR

mmind = where(m09dist NE 0., nmm) & print,nmm

cgOplot,duse[0,mmind],m09dist[mmind],psym=14,color=nirexcolor,$
        symsize=0.35*1.5
m09bad = where( abs(duse[0,mmind]-m09dist[mmind]) GT 3.5, nmbad)

print,'NMBad: ',nmbad
;; cgOplot,duse[mmind[m09bad],0],m09dist[mmind[m09bad]],psym=14,$
;;       color='forest green',symsize=0.75*symper

print,emaf[mmind[m09bad]].cnum
print,emaf[mmind[m09bad]].c
print,m09dist[mmind[m09bad]]


print,'NIREX Offset: ',m4_stat(-duse[0,mmind]+m09dist[mmind])

;;============================================================================
;; Plot masers 
restore,'./emaf_paper/training_set.sav',/ver
nt = n_elements(training)

  mas_ind = lonarr(nwc)
  mas_dist = fltarr(nwc)

FOR kk=0,nt-1 DO BEGIN
   ;; Select by distance
   gcirc,2,training[kk].l,training[kk].b,emaf.l,emaf.b,dis
   delv = abs(emaf.vlsr-training[kk].vlsr)
   
   ;; ind = WHERE(dis LE 0.25*3600. AND delv LE 10., nrt)
   ind = WHERE(dis LE 0.25*3600. AND delv LE 10., nrt)
   print,training[kk].name, nrt
   
   IF nrt NE 0 THEN BEGIN
      
      FOR jj=0,nrt-1 DO BEGIN
         
         IF abs(duse[0,ind[jj]]-training[kk].d) GT 1.02 THEN BEGIN
            print,emaf[ind[jj]].cnum,duse[0,ind[jj]],training[kk].d
            col = masercolor
            oploterror,duse[0,ind[jj]],training[kk].d,duse[1,ind[jj]],$
                       training[kk].dm,/lobar,psym=3,errcolor=col,$
                       color=col,symsize=0.75*symper
            oploterror,duse[0,ind[jj]],training[kk].d,duse[2,ind[jj]],$
                       training[kk].dp,/hibar,psym=3,errcolor=col,$
                       color=col,symsize=0.75*symper
         ENDIF
         cgOplot,duse[0,ind[jj]],training[kk].d,psym=17,$
                 color=masercolor,symsize=0.35*1.5
      ENDFOR
   ENDIF
ENDFOR

al_legend,['Maser','GRS','NIREX'],box=1,/top,/right,/clear,$
          psym=[17,16,14],symsize=[1.5,1,1.5]*.35,charsize=0.9*cs,$
          color=[masercolor,'Opposite',nirexcolor]


myps,/done
myps,'./emaf_paper/plots/distance_comparison_paper_f13b_bw.eps',$
     xsize=3.1,ysize=3.1,/bw


plot,[0],[0],charsize=cs,/nodata,xr=range,yr=range,/isotropic,/xst,$
     xtit='Derived Heliocentric Distance  [kpc]',$
     xmargin=[8,1],ymargin=[4,1],yst=5

cgAxis,yaxis=0,yr=[-1,1],ytit='d!dComp!n - d!dDPDF!n  [kpc]',/save,charsize=cs
cgAxis,yaxis=1,ytickformat='blank_axis'

duse = full_duse
emaf = full_emaf
duse = duse[*,gind]
emaf = emaf[gind]

cgOplot,duse[0,*],emafgrs.dist-duse[0,*],psym=16,symsize=0.35

duse = full_duse
emaf = full_emaf

cgOplot,duse[0,mmind],m09dist[mmind]-duse[0,mmind],psym=14,color=nirexcolor,$
        symsize=.35*1.5

FOR kk=0,nt-1 DO BEGIN
   ;; Select by distance
   gcirc,2,training[kk].l,training[kk].b,emaf.l,emaf.b,dis
   delv = abs(emaf.vlsr-training[kk].vlsr)
   
   ;; ind = WHERE(dis LE 0.25*3600. AND delv LE 10., nrt)
   ind = WHERE(dis LE 0.25*3600. AND delv LE 10., nrt)
   print,training[kk].name, nrt
   
   IF nrt NE 0 THEN BEGIN
      
      FOR jj=0,nrt-1 DO BEGIN
         
         IF abs(duse[0,ind[jj]]-training[kk].d) GT 1.02 THEN BEGIN
            print,emaf[ind[jj]].cnum,duse[0,ind[jj]],training[kk].d
            col = masercolor
            oploterror,duse[0,ind[jj]],training[kk].d,duse[1,ind[jj]],$
                       training[kk].dm,/lobar,psym=3,errcolor=col,$
                       color=col,symsize=0.75*symper
            oploterror,duse[0,ind[jj]],training[kk].d,duse[2,ind[jj]],$
                       training[kk].dp,/hibar,psym=3,errcolor=col,$
                       color=col,symsize=0.75*symper
         ENDIF
         cgOplot,duse[0,ind[jj]],training[kk].d-duse[0,ind[jj]],psym=17,$
                 color=masercolor,symsize=.35*1.5
      ENDFOR
   ENDIF
ENDFOR


al_legend,['Maser','GRS','NIREX'],box=1,/top,/right,/clear,$
          psym=[17,16,14],symsize=[1.5,1,1.5]*0.35,charsize=0.9*cs,$
          color=[masercolor,'Opposite',nirexcolor]

myps,/done


;; plothist,-duse[mmind,0]+m09dist[mmind],bin=0.5,xtit='M09 - DPDF [kpc]',$
;;          ytit='N per 0.5-kpc bin'

e = emaf[mmind]
d = duse[0,mmind]
m = m09dist[mmind]
n = n_elements(mmind)

cgPlot,e.l,-d+m,xr=[50,0],$
       xtit='GLON [deg]',ytit='M09 - DPDF [kpc]',psym=16,/nodata
vline,/h,1.0,color='blk4',linestyle=2,thick=3
vline,/h,-1.0,color='blk4',linestyle=2,thick=3

for kk=0,n-1 do begin
   
   cgOplot,e[kk].l,-d[kk]+m[kk],psym=16,color=d[kk] LE 4 ? 'green' : 'black'
   
endfor

END
